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We use the perturbation method to approximately solve subdiffusion- 
reaction equations. Within this method we obtain the solutions of the 
zcroth and the first order. The comparison our analytical solutions with 
the numerical results shown that the preturbation method can be useful to 
find approximate solutions of nonlinear subdiffusion-reaction equations. 
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1. Introduction 

There are a lot of differential equations which have solutions only in a 
very few special cases and their general solutions remain unknown. Such an 
example is nonlinear differential equations with a fractional time derivative 
which describes the subdiffusion-reaction symmetrical system with two ini- 
tially separated diffusing particles of speices A and B reacting according to 
the formula A + B -> 0(meri) [H El E] 

d a A(x,t) d 2 A(x,t) 

d \ a ' = D d y > -kA(x,t)B(x,t), (1) 

d a B(x,t) d 2 B(x,t) 

d \ a = D d y - kA(x, t)B(x, t) , (2) 



* Presented at 24th Marian Smoluchowski Symposium on Statistical Physics 



(1) 



2 



LEWANDOWSKA PRINTED ON JANUARY 18, 2012 



where A and B denote the concentrations of diffusing particles of species 
A and B, respectively, D — their subdiffusion coefficients (the same for 
both substances), k is the reaction rate constant and the Caputo fractional 
derivative d a f(x)/dt a is defined as [I] 

d a f(t) _ 1 ft., f {n \t>) 



dt a T(n-a)Jo (t-t') 1+a - n 

f( n ) denotes the derivative of natural order n and n— 1 < a < n. We assume 
that the substances are separated from ech other at an initial moment and 
the initial conditions are 

Itl: B ^>-{°c;„ itl: «> 

The symmetry of the system gives A(—x,t) = B(x,t). 

Since the general method of solving fractional subdiffusion-reaction equa- 
tions has not been found yet, one usually uses various approximations, such 
as the quasistationary approximation [3] or the scaling method [2j. In this 
paper we present an idea to solve differential equations by means of the per- 
turbation method. The perturbation method is usually used in nonlinear 
equations in the following form 

g = T(f) + eQ(f) , (5) 

where e is a dimensionless small parameter and it is assumed that the solu- 
tion of the equation 

f=^(A), (6) 

is known. Then, the solution of Eq. ([5]) is a power series with respect to the 
parameter e 

oo 

/ = E * n fn ■ (7) 
n=0 

However, the perturbation method cannot be used directly in Eqs. (pQ) 
and ([2]). The reason is that the system described by Eqs. ([T]), ([2]) and dU 
has its own structure where the depletion zone, the reaction region and 
diffusion region occur. As is shown in [3] within the reaction region the 
terms on the right hand side of Eq. ([T|) (or ([2])) are comparable to each 
other and do not fulfil the assumptions of the perturbation method. Let 
us note that the reaction term kAB and the diffusion term Dd 2 {A, B}/dx 2 
have different scaling properties [2]. We use this scaling property and we 
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change the variables {x and t) to the dimensionless variables in such a way 
that the reaction term will be small compared to the diffusion term. Then, 
we need to solve the transformed equations by means of the perturbation 
method, and next we return to the original variables. We expect that the 
perturbation method works for the transformed dimensionless equations for 
e<Cl but it is not obvious that the dimensional approximate solutions are 
accurate. To check the correctness of these solutions we will compare them 
with numerical solutions to Eqs. ([T]) and In our paper we will find the 
solution of the zeroth and the first order of Eqs. ([T]) and As far as we 
know the perturbation method has not been yet applied to solve fractional 
subdiffusion-reaction equations. We add that the preturbation method was 
appiled to solve normal diffusion-reaction equations with one static and one 
mobile reactants [5]. 



2. The perturbation method 

Firstly, we transform the equations (JTJ) and ((2|) to the dimensionless form 
using the substitutions 



X — pX g , t — Tt s , 



(8) 



where p and r denote the dimensionless position and time, x s and t s are 
constants of the dimension of space and time, respectively. We obtain 



where 



a(p,r 



d a a(p,r) d 2 a(p,r) 

d a b(p,r) _ d 2 b{p,r) 
dr a ~ dp 2 

A(px s ,Tt s ) 



■ ea(p,r)b(p,T) , 
ea(p,T)b(p,r) , 



Co ' 
x s = ^Dt% 
with the initial conditions 

■<*">={ S; :i° ; 

and the boundary ones 

o(— oo,t) = 6(oo,r) = 1 , 



b( P ,r) 



B(px s ,Tt s ) 



Co 



e = kt a s Co , 



b(p,0) 



0, p<0, 

1, p>0. 



da(p,r) 



dp 



db(p,T) 



p=0 



dp 



(9) 
(10) 



(11) 



(12) 



(13) 



p=0 
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We assume that the solutions to Eqs. ([9|) and (|10p are given in the form 
of © 

oo oo 

a(p, r) = £ e n a n (p, r) , b(p, r) = £ e n b n (p, r) . (14) 

n=0 n=0 

Substituting Eqs. (fT4"|) to Eqs. Q and (fTTJ|) and comparing the functions 
of the same order with respect to the variable e occurring on the both sides 
of these equations we obtain the following equations for the functions of the 
zeroth order 



d a a (p,T) _ d 2 a (p,T) 

dr a ~ dp 2 
d a b (p,T) d 2 b ( P ,T) 



(15) 
(16) 



0r a dp 2 ' 

with the initial conditions 

a ( P ,o) = e(- P ) , b (p,o) = e(p) , (17) 



where G is the Heaviside function, and the boundary conditions 

ao(-oo,r) = 1 , f 6 (-oo,r) = 0, 

a (+oo,r) = 0, \ 6 (+oo,r) = l. 



(18) 



For n = 1, 2, 3, ... we obtain the equations for the functions of n-th order 
d a a n (p,T) d 2 a n (p,r) 



dr a dp 2 
d a b n (p,r) d 2 b n (p,T) 



dr a dp 2 



RnM, (19) 
Rn(p,r), (20) 



where 



n-l 

Rn(p,r) = ^a k (p,T)b n _ k {p,T) , (21) 

with the initial condition 

a n {p,0) = b n (p,T) = , n>l, (22) 

and the bonduary ones 

a n (-oo, t) = a n (+oo, r) = , & n (-oo, r) = 6„(+oo, r) = . (23) 

We solve Eqs. (fT5|h ^6|) . (fl9|) and ([20]) by means of the Laplace transform 
method The zeroth order equations (|15p and (TTU|) appear to be subdiffu- 
sion equations without chemical reactions. Equations (|19|) and (|20p become 
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even more difficult to solve when the order of the perturbation method in- 
creases. We will find the exact solutions for a,{ and bi where i = 0, 1. When 
the parameter e is small (e <C 1), the approximate solutions to Eqs. ([9]) 
and (PTOl) are assumed to be 



a = ao + eai , 



b + eb 1 . 



(24) 



3. Approximate solution to subdiffusion reaction equations 

The calulations are arduous, thus the details of these calculations will 
be presented elsewhere. In term of dimensional variables the solutions are 



A ± (x,t) = Af(x,t) + Af(x,t) , 
B ± (x,t) = A ± (-x,t) , 



where 



A (x,t) 
A+ (x,t) 



Co 



2 ^ rn 



i 



i=0 



2 f-^rfi 



ak/2)k\ Vv^DW 2 

x 



i=0 



ak/2)kl 



Dt a / 2 



(25) 
(26) 

(27) 
(28) 



The function Af(x, t) is controlled both by chemical reactions and diffusion. 
It can be treated as a corrections of A^(x,t) which is controlled by the 
diffusion process only. Thus, Af makes samller the particles concentration 
described by the pure diffusion solution and has a proper physical meaning 
for negative values only, and it is expressed by the formula 



Af(x,t) = C 2 kt a F ± (x,t) 



(29) 



when F ± (x,t) < 0, but we put Af (x,t) = when F ± (x,t) > 0, where 

F\ X ,t) -\± Qi (^\ + fpt f-T^V , (30) 

i=0 \ / i=l 



Dt a / 2 



Qi 



l 



(2i)!r(l - (i-l)a) 
1 

i!r(l -a(i/2 - 1)) 



L(i+l)/2j-l 



E 

m=0 



vn\e m , 



(31) 
(32) 
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[u\ denotes here an integer part of a number u ([u\ < u), e$ = 1/4, ef = 0, 
and for m > 2 we get 




-6 -4 -2 2 4 



X 

Fig. 1. Concentration profiles calculated for a — 0.7, D = 0.05, Co = 1, k = 0.001 
and for the times given in the legend; symbols represent the numerical solutions, 
lines — the approximate analytical solutions (|26[) and the zeroth order solutions 
(i.e. for the system without chemical reactions). 

Let us note, that in the functions (|25p - (|34p the parameter e is absent. 
This absence is caused by making the reciprocal inverse transformation of 
the variables (from (x,t) to (p, r) and vice versa). As we have mentioned 
above, we expect that the perturbation method gives the accurate solutions 
for dimensionless equations (|19j) and (j20|) when e <C 1. However, it is 
not clear whether after transforming the dimensionless solutions into the 
dimensional ones these solutions will be correct. In order to ensure this, 
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we compare the solutions obtained by the perturbation method (j26|) with 
the numerical solutions to the sub diffusion-reaction equations (P) and ([2]). 
The numerical procedure of solvnig the subdiffusion equations is presented 
in [31 [7]. In Figs. [T]and [2] we present the comparision between the analytical 
and numerical solutions. We observe that this agreement is reasonably 
accurate for the times presented in the figures. However, as we can see 
in Figs. [T] and [2] the similarity between the approximate solutions and the 
numerical ones decreases with ascending time. The zeroth order solutions 
corresponding to the solution to pure subdiffusion equations do not match 
the numerical solutions accurately but after adding the correction of the 
first order this match is considerably improved. 




-10 -8 -6 -4 -2 2 4 6 



x 

Fig. 2. The same situation as in Fig. Q] but for a — 0.7, D = 0.05, Co = 1, 
k = 0.0001. 



4. Final remarks 

In this paper we find the approximate solution to sub diffusion-reaction 
equations ([T|) and ([2|) by means of the perturbation method. We calculate 
this solution for the zeroth and the first order of the perturbation method 
alone. We also compare our analytical solutions with the numerical ones. 
This comparison is presented in Figs. [T] and [2J The considerations presented 
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in our paper show that the perturbation method can be useful in solving 
nonlinear sub diffusion-reaction differential equations with fractional time 
derivatives. 
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